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We apply theoretically open-loop quantum optimal control techniques to provide methods for the verification 
of various quantum coherent transport mechanisms in natural and artificial light-harvesting complexes under 
realistic experimental constraints. We demonstrate that optimally shaped laser pulses allow to faithfully prepare 
the photosystem in specified initial states (such as localized excitation or coherent superposition, i.e. propa- 
gating and non-propagating states) and to probe efficiently the dynamics. These results provide a path towards 
the discrimination of the different transport pathways and to the characterization of environmental properties, 
enhancing our understanding of the role that coherent processes may play in biological complexes. 



INTRODUCTION 

Recent experimental work has provided evidence support- 
ing the existence of long-lived electronic coherence during 
excitation energy transfer (EET) in photosynthetic complexes 
0][2]. Subsequent theoretical work has then highlighted the 
importance of an intricate interplay of noise and quantum co- 
herence for the efficiency of excitation energy transfer (EET) 
in light harvesting complexes during photosynthesis and iden- 
tified the crucial building blocks that underly this interplay 
ll3l fT0l . While computer simulations and analytical work al- 
low us to identify and verify the importance of these effects 
in theory, the experimental verification of their relevance in 
actual bio-molecular systems is still outstanding. One path 
forward towards this goal is the application of optimal control 
methods to the bio-molecular quantum dynamics in contact 
with environments with the aim of preparing specific initial 
states and control the subsequent dynamics. This would then 
allow for the generation of dynamical behaviour and hence 
signals that can lead to the biggest discrepancies between al- 
ternative theoretical hypotheses. Indeed, here we employ op- 
timal control to develop different strategies that, when exper- 
imentally tested, will allow to enhance our comprehension of 
coherent processes in biological complexes. 

Quantum coherent control drives a quantum system dynam- 
ics towards a specific goal by exploiting quantum coherence 
and interference effects [11 1. Coherent control techniques 
have been proposed for photochemical and photobiological 
processes - see Refs. lfl2lTT7l for an overview on this topic. 
For instance, shaped light has been used to discriminate spec- 
troscopically indistinguishable biochromophores through se- 
lective fluorescence depletion [18] and also to control energy 
flow in bacterial photosynthesis by increasing the amount of 
a triplet state signal in comparison to a singlet state signal 
in carotenoids fl9l . The first evidence of control of exciton 
states in light-harvesting systems was presented experimen- 
tally in Ref. |20l . Moreover, to the best of our knowledge, 
the first attempt to use open-loop optimization to control the 
exciton dynamics of a light-harversting system was theoret- 
ically provided in Ref. GTl . however the control algorithm 
employed there imposed some limitations to their analysis: in 
particular their optimization was limited to state populations. 



To overcome these limitations, we present and apply a re- 
cently introduced optimization algorithm (CRAB) E2ll . In 
particular, we apply quantum (open loop) optimal quantum 
control theory to the dynamics of the electronic excitations 
in the Fenna-Matthews-Olson (FMO), a biological pigment- 
protein complex involved in the early steps of photosynthe- 
sis in green sulphur bacteria [T| [2] [2314251 . Specifically, by 
using the CRAB approach in the context of the FMO com- 
plex i) we achieve general state preparation: this will allow 
us to prepare specific initial states, especially fast and slow 
propagating states exhibiting constructive or destructive inter- 
ference, ii) we explore experimental constraints and imperfec- 
tions (adapted to the expected experimental setup [26]), iii) we 
optimize the difference in signals for different preparations to 
test theoretical hypotheses, and iv) we also discuss optimized 
probing of the system by optimal control of the readout pulses. 
Finally we provide an outlook discussing various possible ex- 
tensions of our approach. 

THE MODEL 

In this section we present the basic ingredient of our the- 
oretical model and fix the notation. The effective dynamics 
of the FMO complex can be modeled by a 7-qubit Hamilto- 
nian describing the coherent exchange of excitations between 
chromophores or sites, i.e. 

7 

Hfrno = huJ l (J t a 7 + H hv oA°j°t + ^ +Cr r) 

where aj (crj) are the raising (lowering) operators for site 
j, ftuij is the local site excitation energy, and vjj denotes the 
hopping rate of an excitation between the sites j and I - see 
Ref. for more details about this model. In the site basis, 
we follow [27 1 and employ the Hamiltonian matrix elements 
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where the zero of energy has been shifted by 12230 cm -1 
for all sites, corresponding to a wavelength of = 800 nm (all 
numbers are given in units of cm -1 = 1.988865-10 -23 Nm = 
1.2414 10 -4 eV) - see Fig. [T] In a first approximation, and 
with the aim of identifying the main transport paths, the dissi- 
pation and dephasing caused by the surrounding environment 
are modeled by the following local Lindblad terms 
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C dlss { P ) = r,-[-{o-/o7,p} + 2a-pa+] (1) 

3=1 
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Cdeph(p) = Yl][-Wt^j ,P} + 2(J t a jP' J t (T J^ ' (2) 
i'=l 

with Tj and jj being the dissipative and dephasing rates at 
the site j, respectively. In the following, we will consider the 
case in which the dephasing and dissipation rates are equal for 
all sites and labeled. As in Ref. 0, we denote the common 
rates as 7 = jj and T = Tj = 5 x 10~ 4 ps -1 for site j 

1 7. The latter corresponds to the measured lifetime of 

excitons which is of the order of 1 ns. Finally, the transfer 
efficiency into the reaction center is measured in terms of the 
population in the 'sink', numbered 8, which is populated by 
an irreversible decay process (with rate T S i nk ) from the site 3, 
as described by the Lindblad term 

C sink (p) = r sinfe [20-^0-3 papers - {(T^ er 8 a£ ', p}] (3) 

where V S ink ~ 6.3 ps -1 (note that H ~ 5.3 
cm -1 ps). The transfer efficiency is given by p s ink(t) — 
2r sinfe J Q p 33 (t')dt', with p 33 (f) being the population of site 
3 at time t'. In the inset of Fig. [2] we show the behavior of 
Psink as a function of the dephasing rate 7, at time t ~ 10 ps, 
when one excitation is initially injected in site 1. 

In order to describe the coupling between the FMO com- 
plex and a short laser pulse, typically used in the laboratory to 
irradiate it |[Tl l23ll24ll2~7l . we add also a semiclassical time- 
dependent Hamiltonian term, HpMo-iaserit), which in ro- 
tating wave approximation takes the form 

7 

Hi as (t) = -J2 ik ■ eE(t) e~ iuit a+ + h.c. (4) 

i=l 

where fa is the molecular transition dipole moment of the in- 
dividual site i ESI , e and uii are, respectively, the polarization 
and the frequency of the field, and E(t) is the time-dependent 
electric field. In the following, we assume E(t) having the 
form 

E{t) = E f(t) , 

with E = 15 D- 1 cm' 1 - 9 • 10 7 V/m (where in SI 
units the Debye is given by D ~ 3.34 • 10~ 30 C ■ m), and a 
time-dependent modulation 

_ er ' 2 $ 1 + J2k=i A k sin(V fc t) + B k cos(v k t) 

A(t) i + Er=ii^i + i^i 

(5) 
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HG. 1. FMO energy level structure where \i) denotes a single exci- 
tation in site i. The states |±) are the symmetric and anti-symmetric 
superpositions of the states |1) and |2). The green and red bubbles 
identify the fast and slow transport path, as detailed in Ref. |9 1 where 
the interplay between different transport pathways in FMO dynamics 
was discussed at length. The main effect of the inclusion of dephas- 
ing noise is the opening of an incoherent relaxation channel from 
level I — } to level |+) (red wiggled line) and therefore the effective 
suppression of the coherent oscillation between level |— ) and sites 
6-7-4 that dominates the coherent dynamics and is responsible of the 
very slow transport once the sink population has reached 50% (initial 
population held in site |+)). The proposed quantum control strate- 
gies efficiently probe this dynamical model. 

with a ramp factor \{t) = i +5 [ e 200(t-T)/T +e -200t/T] (fjuch 
that /(0) - f(T) - 0), and A k , B k , and v k = 2irkr/T be- 
ing parameters to be optimized by using the method described 
below, where r is a random number, T is the time at which 
we want to prepare, for instance, some desired state. More- 
over, we vary also the angles 8 and <j> of the polarization axis 
e, with respect to the dipole moment of site 1, i.e.6> = 6\ + A9 
and 4> = <fii + A(/>, with 61 and <fii describing the orienta- 
tion of the site-1 dipole moment, and A8 and Ac/) being some 
free parameters. The dipole moments of the 7 FMO chro- 
mophores along the three reference axes are (in Debye D, 
where 1 D - 3.34 • 10~ 30 C ■ to): 
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In the following, we choose to = 7 and T = 250 fs and also 
chose the carrier frequency uii as an additional free parameter 
in the optimization. We have also worked with higher values 
of m, up to to = 25, which resulted in small fidelity enhance- 
ments (a few %) but resulted in longer optimization times. 
As direct single site addressing is not possible in the FMO 
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complex due to the strongly overlapping lines we will apply 
quantum optimal control tools as described more carefully in 
the next section to shape the laser pulse, i.e. E(t), and prepare 
the system in a desired physical state and control its transport 
dynamics. Let us point out already here, that the parametriza- 
tion of the laser pulse has been chosen to include some ad- 
ditional constraints, related to typical constraints that can be 
expected in future experiment in this direction. In particular, 
we limit the laser power both due to experimental reasons and 
in order to avoid damage to the sample or cause strong satura- 
tion. We also impose a limit on the time resolution of the laser 
pulse, by means of a limit to the spectral width of the pulse, 
to ensure that it does not exceed 10 fs, since faster modulation 
of the laser is difficult to control experimentally. 

OPTIMAL CONTROL: BACKGROUND AND METHOD 

Coherent control of exciton states in light-harvesting sys- 
tems was demonstrated experimentally in EOl and theoret- 
ically in BTll . In Ref. [20] feedback-optimized femtosec- 
ond pulses are applied to the LH2 antenna complex from 
Rhodopseudomonas acidopila and to a bioinspired artificial 
dyad molecule, in order to control the efficiency of the light- 
harvesting dynamics. Specifically, they optimize the branch- 
ing ratio of energy transfer between intra- and intermolecular 
channels in the complex's donor-acceptor system and obtain 
an enhancement of about 30% in the LH2 system and about 
10% in the artificial dyad molecule, by shaping the pulses em- 
ploying feedback in an iterative learning loop scheme. 

In all mentioned experimental demonstration of control 
of photochemical and photobiological processes lfT2HT7l [20l . 
closed-loop optimization by evolutionary algorithm was ap- 
plied l29l . The procedure consists of three basic components: 

1) a pulse shaper, generating the pulse shape to be tested, 

2) the experiment, generating the feedback signal by pump- 
probe spectroscopy, 3) a computer, running the learning al- 
gorithm and driving the optimization. Hence, the closed-loop 
optimization proceeds along the following steps: i) a random 
guess of a set of pulses is shaped through the pulse shaper 
and then tested on the sample; ii) then the feedback signal is 
evaluated and used to start an evolutionary genetic algorithm 
(based on selection of 'parents', 'mutations', 'recombination', 
and 'generation' of new sets of pulse shapes); iii) a new set of 
pulses is obtained through the pulse shaper and applied to the 
sample. These steps will be repeated until the optimization 
has converged by following the so-called learning curve. 

Despite the interesting applications of this technique, 
closed-loop optimization tends to be effective only for pop- 
ulation control, while coherent control experiments cannot be 
performed because of the inherent shortcoming of transient 
absorption (TA) spectroscopy. One way to overcome this is- 
sue could be to use 2D electronic spectroscopy in order to get 
information about the phase, which is necessary as a feed- 
back signal in coherent control experiments. Here, we solve 
this issue, using open-loop control approach: one first numeri- 



cally optimizes laser pulses via numerical simulations and the 
applies them to the sample obtaining the desired result, for 
example, the experimental preparation of the sample in some 
desired state with very high fidelity. 

The main advantages of the open-loop approach with re- 
spect to the closed-loop approach are two-folds. In the lat- 
ter, the pulses are often very complex, highly structured, and 
very demanding to interpret: it is usually rather difficult to 
understand the real physical effect of such series of consecu- 
tive pulses on the system. This limits the understanding of the 
physical processes underlying a certain biological behaviour. 
Moreover, repeated closed-loop experiment rarely result in the 
same genetic algorithm-driven learning curve, increasing the 
difficulties of the analysis of the final optimal series and of the 
error estimation. On the other hand, in the open-loop scheme 
the optimally-shaped pulse is well determined and can be ap- 
plied on the sample repeatedly to increase the signal to noise 
level and to compare the output feedback when changing the 
applied phase. It is then easier to find the explanation for the 
response of the system and the physical mechanism behind 
it. This problem becomes even simpler when the CRAB op- 
timization is used as, as explained below, it results in optimal 
but very simple, robust and structured pulses. Let us stress 
however that the open loop technique can be applied directly 
only when the system parameters (e.g., Hamiltonian and envi- 
ronmental noise) are sufficiently well known and this, indeed, 
makes the closed-loop approach more feasible for several bi- 
ological systems in which this information is not yet acces- 
sible. Here recently developed methods for quantum process 
tomography applied to multi-chromophoric systems, provid- 
ing the decoherence of the system, the density matrix, as well 
as the Hamiltonian parameters [ 30 1 will be able to assist our 
open-loop approach. This is even more so as the robustness of 
the pulse shapes that we obtain from our open-loop technique 
allow us to use open loop scheme even if the system details are 
not well measured. Quite reasonably, the best scheme would 
be a combination of these two approaches. Indeed, one could 
use, for instance, closed-loop control as a means to obtain in- 
formation about system parameters, e.g. if one optimizes the 
pulse in the experiment then one can chose this pulse and find 
out theoretically for which system parameters it reproduces 
the experimental data. Then, varying the target that is to be 
optimized, we obtain different pulses and can repeat this in- 
vestigation; each time one obtains useful information about 
the system parameters. Working out such a programme might 
be useful as it might allow an alternative to tomography. Once 
we apply this procedure, by using the obtained information 
about the system, an open loop control can be then success- 
fully applied. 

The first work where open-loop optimal control in FMO 
complexes is reported in |2TI . The authors investigate the con- 
trol of exciton dynamics in FMO complex, by using polarized- 
shaped pulses optimized by means of a derivative functional 
equation for the target function. In particular, they use shaped 
pulses to optimize the excitation energy localization in a sin- 
gle chromophore of the FMO complex (site 7). 
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Here we extend and improve this result using the recently 
developed CRAB optimization technique (described below) 
targeting different initial states, to investigate different trans- 
port and decoherence processes in FMO. Moreover, we are 
able to consider faster processes that are more robust against 
decoherence as we optimize pulses of a few hundred fem- 
tosecond length. This allows us also to perform coherent con- 
trol, i.e. preparation of a coherent superposition state, as well 
excitation energy localization. Concerning the case of sin- 
gle site preparation, it is not feasbile to perform a clear com- 
parison with the results of Ref. ED . since they studied the 
FMO antenna complex present in the Chlorobium Tedidum 
bacterium, while here we investigate the FMO complex found 
in a slightly different bacterium called Prosthecochloris aes- 
tuarii. However, from a general perspectives, it seems that our 
approach allows us to get much higher fidelity, and, more im- 
portantly, to use sensibly shorter laser pulse lengths (250 fs, 
while their optimal pulse are 600 fs long). These difference 
appear to be even more relevant and crucial when one wants 
to prepare a coherent state in presence of strong dephasing 
noise. 

It should be noted that, since we are considering short pulse 
durations, neglecting the double-exciton states is a good ap- 
proximation, also according to the theoretical results in Ref. 
lETl . Besides, the efficiency of the CRAB algorithm allow us 
to investigate more deeply the issue of random orientations of 
the FMOs in the sample by considering very large samples 
(10 4 FMOs), while in Ref. ETI . due to computational limita- 
tions, only some preliminary results were reported. In fact, we 
will demonstrate that even a partial orientation of the samples 
by means of an external field combined with optimal pulses 
will improve experimental results significantly. Finally, we 
show also how to optimize the probe of the system to improve 
the experimental results even more. 

To achieve all this here, we use the Chopped RAndom 
Basis (CRAB) optimization, introduced in Ref. Il22l . to op- 
timize a specific figure of merit, e.g. population in some 
localized or delocalized state at some final time or the fi- 
nal fidelity with respect to a target state •F(T), by vary- 
ing the control field entering into the Hamiltonian term in 
Eq. (ffl. Introducing the control field parametrization given in 
Eq. Q the functional becomes a multivariable cost-function 
J-(A8, A<f>, Ldi, to, a, {Ak}, {Bk}) on which any standard 
minimization method can be applied. We started with O(10 3 ) 
different initial random configurations and applied a direct 
search algorithm, which does not compute gradient nor Hes- 
sian, to find the function minimum OTI . To minimize an ill- 
variable function, the Nelder-Mead algorithm starts defining 
a M + 1 dimensional polytope and then, in its simplest im- 
plementation, moves it replacing the worst point with a point 
reflected through the barycenter of the other M points, result- 
ing in a (local) minimization of the function. We used the 
Subplex variant of the Nelder-Mead algorithm, which applies 
the same algorithm to different subspaces to improve the con- 
vergence l32l . The CRAB optimization strategy introduced 
above allows -as we shall show in the following- to find the 
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FIG. 2. e a versus dephasing rate 7 (in units of ps~ ), for different 
prepared initial states. Specifically, our goal is to prepare the states 
B and D (i.e., a — B, D). As described in the text, these quantities 
are given by e s = 1 - {+\p\+} = 1 - P11 + P22 - »[pi 2 ] and 
Ed = 1 — ps,5 — pe,6 — pr,r- Inset: Transfer efficiency vs. dephasing 
rate 7 (uniform for all the sites) at t — 10 ps, when one excitation is 
initially in site 1, 

optimal pulses to extremize the desired figure of merit: it is 
efficient and versatile as it does not need any analytical solu- 
tions of the system dynamics, it does not compute gradients 
and can be easily adapted to different figure of merits. More 
importantly, it includes already many experimental constraints 
such as the finite band-width and power of the control pulses. 
Finally, we mention that optimal control pulses are quite ro- 
bust with respect to system parameters perturbation and noise 
up to a few percent, as shown in the literature in different sce- 
narios (see e.g. [33, 34]). This robustness arise from the fact 
that the optimal dynamics lies in a minima of the functional to 
be extremized, thus first order perturbations vanish. 

INITIAL STATE PREPARATION 

In Ref. J51, it has been found that the coherent transfer 
of the electronic excitation energy in the FMO complex takes 
place essentially through two different pathways: one medi- 
ated by the state |+), which is shifted towards site 3 and leads 
to very fast transfer to the reaction center, and a second one, 
involving the sites 5, 6 and 7, which is comparatively slow 
because of the energy gap with the site 3 and because the 
excitation suffers many coherent oscillations between those 
sites before reaching site 3. Indeed, the presence of dephas- 
ing noise assists the transport because, on the one hand, it 
opens up a new additional pathway, i.e. incoherent tunneling 
between the state |— ) and |+), and, on the other hand, it par- 
tially suppresses the transition from |— ) to sites 5, 6, 7, and 
also leads to fast incoherent oscillations between those three 
sites before reaching the sink. Motivated by these results, 
we apply the CRAB optimization to find the optimal pulses to 
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FIG. 4. Modulus of the ensemble average of the elements of the FMO 
density matrix in the site basis (over 10 4 samples), after applying the 
optimal laser pulse to prepare the system in the state B. We introduce 
also static disorder, set to 1% (left) and 100% (right), and we get a 
fidelity for the desired state equal to 89% and 49.8%, respectively. 



FIG. 3. Probability distribution of the quantity eb for the preparation 
of the state B, in the presence of dephasing with rate 7 ~ 1 ps , 
by using the optimal and a standard Gaussian laser pulse. We con- 
sidered a sample of 10 4 FMO complexes to get enough statistics. 
Inset: shape of the optimal pulse, optimized to prepare the state B 
in 250 fs. Notice that the laser pulse amplitude is always shown in 
units of D" 1 cm" 1 ~ 6 • 10 6 V/m. Finally, the optimal val- 
ues of the other free parameters are AO = 2.5, A<fi = 7.6, and 
uii — 121.76 cm~ , and we find eb = 0.20, without averaging. 
Interestingly enough, by considering a simpler linear shape (dashed 
line), we find a quite similar error, i.e. eb = 0.24, by showing the 
'robustness' of the pulse shape to prepare the state B. 

selectively prepare a state \D) with maximum probability of 
finding the electronic excitation in the sites |5), |6), |7) (dark 
or non-propagating state) and \B) = |+) = (|1) + \2))/\/2 
(bright or propagating state). We consider the full model for 
the FMO complex (described before), but in absence of reac- 
tion center, since it is usually the case in the current experi- 
ments on this light-harvesting system. However, the presence 
of the sink would not affect the state preparation process since 
the laser pulse is applied for a very short laser pulse. This is 
also confirmed by the analysis on the transport pathways dis- 
cussed below. In particular, we apply the CRAB strategy to 
prepare the FMO complex (initially in the ground state) in the 
desired state after t = 250 fs during which we excite the sys- 
tem with a laser pulse. We maximize the probability of finding 
the excitation in the desired state (B or D), that is our figure 
of merit, by minimizing the following quantities 

e B = l- (+\p\+) = 1 - Pn + p22 - jft[p 12 ] , (6) 

£ D = 1 - p 5i 5 - p 6 ,6 - P7,7 • (7) 

In Fig. [2] results of the optimization are shown as a function 
of the dephasing rate present in the system. As it can be seen, 
increasing the amount of dephasing in the dynamics, the error 
e a increases in both cases. 

As remarked earlier, the so achieved pulses are quite robust 
against changes in system parameters such as the environmen- 
tal noise level. To demonstrate this we apply the following 



procedure. We find the optimal laser pulse in the presence of 
dephasing with rate 7 ~ 1 ps -1 . Then, we repeat the opti- 
mization with different values of dephasing and find, in gen- 
eral, only slightly different pulse shapes that result in values 
of error e a , which differ only slightly, less than 10 -2 , from 
those obtained for the optimal shapes we found for dephasing 
7 ~ 1 ps" 1 . In other words, the same optimal pulse can be 
used to prepare some desired state irrespective of the strength 
of the dephasing noise in the dynamics is. These pulses are 
shown in the insets of Figs. [3] and [5] Let us also point out that, 
by considering a simpler linear pulse shape (interpolating the 
optimal pulse, see dashed line in the inset of Fig. [3}, one ob- 
tains a very similar error, whose difference is less than 0.05 in 
both cases (D and S states). As the optimal pulse shapes de- 
pend only marginally on the noise level, the results in Fig. 
12 can become a tool to obtain information about the noise 
strength in the FMO dynamics by measuring the quantity e a . 
Notice also that, removing the experimental constraints on the 
pulse shape - namely the limit on the maximal pulse intensity 
and the time resolution - it is possible to find near-unit fidelity 
for any state and any value of the dephasing rate. For instance, 
if one allows twice as high field strength, an improvement of 
even several % for the fidelity can be obtained. 

In realistic samples, the orientation of the FMO complexes 
does exhibit significant disorder. Hence it is crucial to con- 
sider the effect of the effect of random orientation of the FMO 
complexes which leads to an ensemble average over the ran- 
dom distributions. In particular, we consider two extreme 
cases in which we add 1% and 100% of random disorder, i.e. 
= op t + VSi and cj) = <f> opt + r)s 2 , where 8 opt and <f> opt 
are the optimal values, Si 2 are random numbers in the range 
[0, 27r] and r\ = 0.01,1 respectively. As illustrated by the 
results in Figs. [3] and [5] the corresponding distributions are 
fairly distinguishable. Indeed, for the state B analyzed in Fig. 
[3] the ensemble averages of eb are 0.207 and 0.751, respec- 
tively for 1% and 100% of orientation disorder, in the case 
of the optimal pulse, while they are 0.793 and 0.904 in the 
Gaussian one, respectively. On the other hand, for the state 
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FIG. 5. Probability distribution of the quantity Ed for the preparation 
of the state D, in the presence of dephasing with rate 7 ~ 1 ps _ , 
by using the optimal and a standard Gaussian laser pulse. As be- 
fore, since in the lab one has a sample of FMO complexes in random 
orientations, we plot these pdfs for two extreme cases in which we 
add 1% and 100% of random disorder to the two angles defining the 
orientation of the FMO complex. Inset: shape of the corresponding 
optimal pulse. The optimal values of the other free parameters are 
A6 = 3.09, A0 = 3.76, and uji = 504.46 cm" 1 , and we find 
Ed — 0.09, without averaging. By considering a much simpler con- 
stant laser field (dashed line), one still gets a very good preparation 
of the state D, i.e. s D = 0.10. 



D analyzed in Fig. [5] the distribution averages are, respec- 
tively, 0.091 and 0.531 in the case of the optimal pulse, and 
0.663 and 0.634 in the Gaussian one. Moreover, the modu- 
lus of the ensemble average of the elements of FMO density 
matrix in the site basis, after the state preparation by means 
of the optimal pulse and the Gaussian one, are shown in all 
the analyzed cases (in presence of static disorder) in Figs. [4] 
[6] [7] Furthermore, we measure the fidelity F(p, a) between 
the desired state a and the one p achieved by optimal control 
as F(p, a) = Ti[y/ yfa p^fa]. Notice that, in the case of state 




FIG. 6. Modulus of the ensemble average of the elements of the FMO 
density matrix in the site basis (over 10 4 samples), after applying a 
standard Gaussian pulse. We consider also the presence of static 
disorder, which is set to 1% (left) and 100% (right), respectively. 




FIG. 7. Modulus of the ensemble average of the elements of the FMO 
density matrix in the site basis (over 10 4 samples), after applying 
the optimal laser pulse to prepare the system in the state D. Again, 
the static disorder is set to 1% (left) and 100% (right) and we get a 
fidelity for the desired state equal to 86% and 63.1%, respectively. 

D, since the quantity is defined only in terms of the total 
population in site 5, 6 and 7, in order to calculate the fidelity, 
we need to specify a particular goal state a. For simplicity, we 
choose a target state with 70% population in site 5, 25% in site 
6, 5% in site 7, and vanishing coherences. This distribution of 
populations is similar to the one obtained when minimizing 
the quantity S£>. In Fig. [4] we notice that, if one will be able, 
in the lab, to orientate very well the FMO complexes in the 
sample, the state B is prepared with very high fidelity. Addi- 
tionally, when the FMO complexes are completely randomly 
oriented in the sample, the fidelity is much lower. Indeed, one 
observes a large amount of population in site 4 which is al- 
most in resonance with the site 1 and 2. However, this does 
not affect the transport pathway discrimination (as shown be- 
low), because the population in site 4 goes quickly to the site 
3 as well as the state B. A similar behaviour is observed for 
the preparation of the state D. Nevertheless, in this case the 
high values of the fidelity are more robust to the introduc- 
tion of orientation disorder, since one wants only to populate 
the high-energy sites 5, 6 and 7 neglecting any control of the 
phase terms. Instead, in Fig. |6j one clearly observes that 
a gaussian laser pulse is not able to selectively prepare spe- 
cific states. Furthermore, when 100% of orientation disorder 
is considered, the gaussian pulse does almost equally popu- 
late all the sites with vanishing cross-term correlations (com- 
pletely mixed state). Finally, it is worth pointing out that an 
experimental partial orientation of the sample will be able to 
sensibly enhance these fidelities, as shown below in the sec- 
tion about the orientation. In the following, we will use these 
optimal laser pulses, found to prepare those selected states, 
in order to investigate the different transport pathways in the 
FMO complex dynamics. 

STATE PREPARATION DEPENDENT TRANSPORT 

We now compare the transfer efficiency into the reaction 
center corresponding to the two different initial states that we 
have introduced before, including the Lindbladian term in Eq. 
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FIG. 8. Transfer efficiency as a function of time (ps) for different 
initial states, i.e. \D), | — } and \B), for two different values of de- 
phasing rate 7 = and 7 ~ 1 ps -1 . Moreover, in one case, we 
consider also an idealized preparation of the state, i.e. e a = 0. In 
the other cases, we consider the more realistic scenario when in the 
lab one prepares similar states applying the optimal laser pulses, in 
absence and presence of dephasing noise. 



Q - see Fig. [8] Notice that now the sink is included in the 
dynamics from the beginning, even if the optimal pulses were 
obtained without sink. Indeed, as discussed also above, it does 
not sensibly affect the state preparation since the pulse is ap- 
plied for a very short time scale compared to the transfer rate 
from the site 3 to the reaction center. 

In the absence of dephasing noise, as expected from the re- 
sults in Ref. [9], the population initially prepared in sites 5, 

6 and 7 is basically trapped and only very slowly gets into 
the reaction center, while the state |+) is almost in resonance 
with the site 3 and the transfer rate to the sink is much faster. 
However, increasing the amount of dephasing, these destruc- 
tive interference effects are reduced, destroying the so-called 
invariant subspaces or dark states (see Ref. |5| for more de- 
tails), and the efficiency discrepancy between the fast and the 
slow pathways decreases, as shown in Fig. [8] Particularly, 
the ratio fast/slow pathway transfer efficiency is about 2.5 for 

7 ~ 1 ps -1 and about 80 for 7 = 0, at time t = 2 ps. These 
results are particularly robust against various possible exper- 
imental inaccuracies: In Fig. [8] we show that, although the 
state preparation by the optimal laser pulse is not perfect, the 
corresponding behaviour is still good enough to distinguish 
the slow and fast transport pathways, compared to the case 
in which the specific states are exactly set in the numerical 
simulation. In other words, it turns out that, even if the state 
preparation is not perfect, one finds a very similar behavior for 
the transfer efficiency. This fact is fundamental to reproduce 
these results in the laboratory since the experimental fidelities 
will be smaller than the theoretical ones. Finally, we also con- 
sider the realistic experimental scenario in which a very large 
ensemble of FMO molecules is studied simultaneously in the 
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FIG. 9. Probability distribution of the transfer efficiency into the sink 
after 2 ps, in the presence of local dephasing with rate 7 ~ 1 ps -1 , 
for the optimal preparations of the state D and the state B, with 
100% of random disorder in the orientation of 10 4 FMO complexes. 
The two distributions are distinguishable with less than 5% of error. 



lab with different random orientations. It turns out that, even 
with 100% of random disorder in the orientations of the FMO 
molecules, the difference between the transfer efficiency af- 
ter 2ps of the 'dark' and 'bright' states is still measurable. 
Indeed, as shown in Fig. [9] the two distributions of transfer 
efficiency are distinguishable with less than 5% of error. The 
ensemble averages are around 0.12 in the case of the state D, 
and 0.39 for the state B. In conclusion, the proposed analysis 
appears to be observable in a real experiment. 



OPTIMAL CONTROL PULSE FOR EXPERIMENTS ON 
DIFFERENTLY ORIENTED SYSTEMS 

In the lab one generally has a sample of many FMO com- 
plexes with random orientations and, hence, the laser pulse 
could be optimized also taking it into account. 

In order to cover almost isotropically all the different orien- 
tations of the photosynthetic system in the sample, we con- 
sider the 20 directions pointing to the 20 vertices of a do- 
decahedron, and we optimize the pulse in order to minimize 
the quantity e a averaged over a sample of 20 differently ori- 
ented FMO complexes. Notice that the computational effi- 
ciency of the CRAB algorithm allows us to consider a much 
higher number of directions and we consider here only 20 just 
for simplicity and for illustration purposes. Indeed, the time 
complexity of such algorithms scales linearly with the num- 
ber of orientations to average over and the optimization in the 
case of 20 directions takes only a time of the order of hours 
on a single standard CPU. 

we show the probability distribution for 



In Fig. 10 



when this new optimal pulse is applied to a sample of 10 4 
FMOs in random orientations. As comparison, we plot also 
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FIG. 10. Probability distribution of the quantity Ed for the prepa- 
ration of the state D, in the presence of dephasing with rate 7 ~ 
1 ps~ 1 , by using the optimal pulses obtained by averaging on a single 
FMO complex with a specific orientation (as in Fig. |5j, by averaging 
on a small sample of 20 (almost isotropic) FMO orientations pointing 
the 20 vertices of a dodecahedron, and on a sample of 21 orientations 
inside a cone with an opening angle of 0. 1 n. Hence, one consider the 
proability distribution of sd when these optimal pulses are applied to 
a sample of 10 4 FMO complexes in completely random orientations 
in the first two cases and, for the last case, randomly oriented inside 
that cone. The shape of the corresponding optimal pulses are similar 
to the one in the inset of Fig. [3] Moreover, we find, for the dodecahe- 
dron case, AO = 2.66, A(j> = 2.61, lji = 277.10 cm" 1 , while, for 
the cone case, A0 = -2.21, A(f> = 1.86, u> x = 274.87 cm' 1 . The 
dashed lines represent the corresponding averaged values, i.e. 0.53 
(single orientation), 0.46 (dodecahedron), 0.29 (cone). The corre- 
sponding distribution widths are also plotted. 



the corresponding probability distribution when we instead 
apply the optimal pulse used in Fig. [5] It turns out that this 
new pulse is somehow more robust and provides an error dis- 
tribution whose width is about one half of the one obtained 
with the pulse optimized with a single FMO system. This 
could suggest a way to improve the state preparation of the 
FMO complex in the real experiment where one has a sample 
containing many FMO molecules in a solution. However, con- 
cerning the preparation of the more delicate state B this tech- 
nique does not provide any noticable improvement. More- 
over, regarding the transport properties, the results are similar 
to those presented in Fig. [9] and the overlap of those two dis- 
tributions is still less than 5% but does not decrease further. 
Notice that this overlap is already rather small to be able to be 
experimentally observed. Finally, we consider an intermedi- 
ate case in which the FMO complexes can be partially oriented 
along a cone-shaped orientations. In particular, we repeat the 
analysis above but for 21 directions inside a cone with a 10% 
opening angle, i.e. O.I71", and we optimize the pulse in order 
to minimize the quantity e a averaged over 21 FMO complex 
evolutions. Then, we apply this optimal pulse to a sample of 
10 4 FMOs in random orientations inside this code and we cal- 



< 1.5 




FIG. 11. Top: Behaviour of the quantity A (defined in the text), in 
units of cm" 1 (~ 0.1 meV) as a function of the angles A8 and 
A(j> defining the orientation of the FMO complex, in the case of 
cui = -1000 cm" 1 and E = 70 D" 1 cm" 1 ~ 42 • 10 7 V/m. 
The difference between the maximum and the minimum value is 
comparable to the thermal energy and the maximum is obtained for 
A6 ~ 1.75 and Acj> ~ 2, which seems to be the preferred orienta- 
tion when the sample is subjected to this constant laser field. Bottom: 
maximum value of the Rabi frequencies fli . eEg versus A6* and A<f>. 
Notice that the values of the detunings are much larger than the Rabi 
frequencies. 



culate the probability distribution for Ed - see Fig. [10] We 
find that the width of the error distribution is further squeezed 
and shifted to smaller errors, i.e. higher fidelities. 



TOWARDS THE ORIENTATION OF THE FMO SAMPLE 

Here, we show how one could try to orientate the FMO sys- 
tems in the experimentally available sample, by using some 
simple classical mechanics arguments. Following Ref. 11351 . 
we model each monomer of the FMO complex as a disk, 
whose mass and radius can be reasonably estimated to be 
equal to M ~ 80 kDa - 15 10" 23 Kg and R = 2 nm (in- 
cluding the protein scaffolding). The corresponding moment 
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FIG. 12. Boltzmann-Gibbs probability distribution in orienting the 
FMOs in the sample, as a function of the angles A8 and A<f>, in 
the case of uji = -1000 cm" 1 and E = 70 D' 1 cm' 1 . The 
dashed area represents the region in which we will sum up all the 
probabilities in the case, for instance, of a 20% cone opening angle, 
as shown in Fig.[T3] Notice, however, that a slightly tilted strip would 
get higher probability for the same width and so further improves 
our results. The Gaussian pulse shape, used as comparison in the 
population behavior, is also shown (dashed line). 



of inertia (with respect to one of its diameters as rotational 
axis) is hence / = \MR 2 ~ 1.125 10" 31 Kg m 2 . The ro- 
tational energy is E = giw 2 , with ui being the angular ve- 
locity in radians per second, i.e. the derivative of the angle 
rotated with respect to time uj = 4f. At room temperature, by 
neglecting the friction due to the presence of a solution and, 
possibly, other more sophisticated effects, just to get an ini- 
tial rough estimation, the time t rot it takes for the system to 
rotate by an angle ir/2 (which is roughly the average angle 
by which a complex has to be rotated) is trivially given by 

6 /is, with E t h 



trot = |J ~ ^ /is, with E t h ~ 25 meV being the ther- 
mal energy. Actually, given that the disk will carry out a sort 
of random walk, the rotation time could be much longer than 
what is estimated here by means of this simple analysis. 

Moreover, we investigate the energy landscape as a function 
of the orientation of the FMO complex, when is subjected to 
a constant electric field Eq polarized along the axis e with 
carrier frequency lui, i.e. the following quantity 



A 



7 i - 
y> \JM_ 

i=l 



UJ; 



UJi 



(8) 



By varying the orientation of the FMO complex in terms of 
the two angles 8 and <fi, the energy difference between the 
maximum and the minimum value is comparable to the ther- 
mal energy and one gets the maximum for A6 ~ 1.75 and 

Acj> 
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2 - see Fig. 
-1000 cm" 1 andE'n 



The other parameters are chosen as 
u>i = -1UUU cm"" ana E = 40 D~ x cm' 1 . Finally, follow- 
ing the simple reasoning above, we explicitly calculate how 
many FMO systems can be oriented around some direction at 



FIG. 13. Sample population ratio of FMOs oriented within a cone 
with some opening angle (20%, i.e. 0.27T, and 50%, i.e. 0.57r), as 
a function of temperature (K), with uji = —1000 cm~ , and Eq = 
70 D" 1 cm , according to the Boltzmann-Gibbs distribution. 



a certain external temperature, since the thermal fluctuations 
will unavoidably try to disorientate the sample. To do this, we 
calculate the probability for the FMO to be oriented at a cer- 
tain angle by assuming that they follow a Boltzmann-Gibbs 
distribution, in the presence of a constant electric field - see 
Fig. 



12 Then, by using the obtained probability distribution, 
we also evaluate how many FMO systems are oriented within 
a cone with a certain opening angle as a function of the tem- 
perature and for different values of opening angle. The corre- 



sponding results are shown in Fig. 13 By decreasing the tem- 
perature, the amount of oriented FMO systems in the sample 
increases and this would give us higher fidelities in the state 
preparation analysis and for the other related results above. 
Therefore, it seems feasible to orientate the FMO complexes 
in the sample to some extent using far detuned laser light. The 




FIG. 14. Modulus of the ensemble average of the elements of the 
FMO density matrix in the site basis (over 10 4 samples), after apply- 
ing the optimal laser pulse to prepare the system in the state B (left) 
and D (right), respectively. We introduce a static disorder, set to 10% 
in both cases, considering that some orientation will be obtained in 
the experiment. The fidelity for the preparation of state B and D are 
74.4% and 85.1%, respectively. 
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typical individual optimal control experiments take of the or- 
der of ~ ps which, importantly, are several orders of magni- 
tude shorter than the time for the FMO complexes to loose 
orientation after the detuned laser beam has been switched 
off. Hence, we can avoid interference between the orientation 
laser and the actual optimal control experiment. Finally, let us 
point out that a partial orientation of the sample will allow one 
to sensibly improve the state preparation results shown above 
see Fig. 14 Indeed, a net improvement is observed, when 



I Gaussian pulse 
I Optimal pulse (21 FMO orientations) 



compared to the 100% disorder case shown in Figs. [4]and[7] 



OPTIMAL PROBE 

In the context of controlling a molecular systems by optimal 
femtosecond laser pulses, a similar approach can be reason- 
ably used to optimize the probe pulse absorption in a pump- 
probe scheme. A successful demonstration of the optimiza- 
tion of the absorbance of the probe pulse by optimal control 
techniques, but based on a derivative functional equation, was 
shown for a prototypical molecular three level system in Ref. 
Il36l . Here, we repeat the analysis above for a probe laser 
pulse applied to the FMO complex by using the CRAB algo- 
rithm. In order to compare our theoretical predictions, e.g. 
in Fig. [9] with the experimental data, since there is no reac- 
tion center in the FMO complex sample used in the lab, one 
has to measure the population in site 3 as a function of time 
and then calculate the corresponding transfer efficiency, as de- 
fined in the model in Eq. p). To do that, a probe laser pulse is 
applied to the sample and the corresponding absorption inten- 
sity is detected. Usually, the probe pulse is a Gaussian laser 
pulse on resonance with the site whose population one wants 
to measure, e.g. site 3. Here, we apply the optimal control 
tools to analyze whether a shaped pulse can detect the site 
population, particularly in site 3, with an higher 'resolution', 
as compared to a simpler Gaussian pulse. Actually, it will turn 
out that, when one considers only a single FMO complex, very 
high fidelity (99%) are already obtained by a gaussian pulse 
oriented along some optimal polarization axis and the pulse 
shaping will not give significant improvements. On the other 
side, if one has a sample of fully randomly oriented systems, 
both gaussian and an optimally shaped pulse bring to very low 
fidelities since the orientation disorder is too strong. However, 
if a partial orientation will be feasible from the experimental 
point of view, the control on the pulse shape sensibly increase 
the probability of successfully probing site 3. Specifically, we 
use the model above for a single FMO complex and we con- 
sider the case where initially all the population is in site 3 and 
we want to find the pulse which is able to detect this popula- 
tion by absorption, i.e. by removing population from the site 
3. To apply the optimal control approach as defined above, we 
use the following error function 




with P33 being the site-3 population. 



(9) 



FIG. 15. Probability distribution of the quantity ep, when the FMO is 
initially prepared with all the population in site 3 and then subjected a 
laser probe pulse for a time interval of 125 fs to detect the site-3 pop- 
ulation, in the presence of dephasing with rate 7 ~ 1 ps - . In par- 
ticular, we consider the case of a Gaussian pulse on resonance with 
the site 3, i.e. uii — cm" 1 , and polarized along the optimal (wrt a 
single FMO) orientation, i.e. AO = 2.92, A<j> = 1.85, and the case 
of an optimal pulse obtained by averaging the error function e p over 
21 directions inside a cone with an opening angle of 0.1 n, whose 
optimal parameters are A9 — 5.61, A(f> — 1.39, uji = 12.63 cm' 1 . 
To get the probability distribution, we apply both pulses to a sample 
of 10 4 samples randomly oriented inside that cone. Inset: shape of 
both (Gaussian and optimal) probe pulses. 



Moreover, the pulse is applied for a time interval of t = 
125 fs. We find that, if there is just one FMO complex in the 
sample, the optimization will simply orientate the probe pulse 
polarization axis along some optimal direction (related to the 
site 3 dipole moment). Apart from the polarization control, 
the further optimization of the pulse shape does not bring fur- 
ther significant gains. Hence a Gaussian pulse will reliably 
detect the site-3 population (absorption efficiency, i.e. 1 — ep, 
of 99%). However, if one applied this Gaussian pulse polar- 
ized along this optimized direction to a sample of randomly 
oriented FMOs, the absorption efficiency shows an almost flat 
probability distribution in the range [0.45, 0.99], which is ac- 
tually a bit higher for smaller efficiencies, i.e. the method 
does not provide us with a more efficient absorption signal 
(compared to the traditional way), and has associated an er- 
ror of about 50%. If we apply the optimal control algorithm 
to find the best probe, averaged over 20 isotropic orientations 
according to the dodecahedron above, this does not improve 
sensibly the results either (data not shown). Finally, we con- 
sider the case in which the sample is partially oriented, partic- 
ularly within a cone of 0.1 7r opening angle, and we optimize 
the pulse along 21 orientations inside this cone, as done above 
for the pump. Hence, we apply this optimally shaped pulse to 
a sample of 10 4 FMOs randomly oriented inside this cone and 
we compare the absorption efficiency to the case of a Gaussian 
pulse, polarized along the optimal orientation obtained for a 
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single FMO. We find that the optimally-shaped pulse gives 
a probability of high absorption efficiency (small e p), which 
is more than twice larger than the one with a Gaussian pulse 
- see Fig. 



15 



Therefore, in presence of structural disorder 
but achieving a partial orientation of the FMOs in the sam- 
ple, both polarization and shape optimizations may enhance 
the probe pulse absorption in a pump-probe scheme, which is 
crucial for an efficient discrimination of dynamical properties, 
like the identification of transport paths discussed here. 



SUMMARY AND PROSPECTS 

In summary, the experimental verification in bio-molecular 
systems of fundamental building blocks of quantum dynam- 
ics in presence of environment [3-10] requires the develop- 
ment of novel experimental tools and theoretical methodol- 
ogy. In this work we have contributed to this effort with the 
demonstration that novel methods from the theory of opti- 
mal control can be combined with ultra-fast laser pulses to 
provide enhanced diagnostic tools suggesting promising new 
routes for experiment. In particular we have introduced and 
applied the CRAB technique for optimal quantum control that 
was originally developed in quantum information science El 
and used it to determine, for realistic experimental parameters 
ll26l . pulse shapes that allow for the preparation of arbitrary 
coherent superpositions with high fidelity. We are also able 
to reduce the impact of the random orientation of FMO com- 
plexes in typical samples and we are able to optimize the read- 
out of the system to maximize state sensitivity. Finally, we in- 
dicated that, in a future experiment applying this approach, a 
comparison of our theoretical calculations with experimental 
data, would allow us to extract information about the largely 
unknown details of the system-environment interaction, thus 
complementing recently proposed methods based on tomog- 
raphy [30]. These methods may support the experimental con- 
firmation that recent models concerning the interplay of trans- 
port processes and environmental noise [3-10] grasp the main 
features of the system dynamics. 

In future work, based on the techniques presented here, we 
are planning to consider more general non-Markovian mod- 
els J37 -40], other bio-molecular complexes and we will ex- 
plore the importance of multiple excitations in the system. 
We will also explore the use of CRAB as a tool for closed- 
loop control in this context. These tools will then be applied 
to provide optimized experimental set-ups that explore ques- 
tions concerning the relevance of quantum coherence and de- 
coherence in the dynamics of bio-molecular systems. This 
includes the quantitative probing of the functional relevance 
of quantum coherence and entanglement, the exploration of 
which would shed further light on the question whether entan- 
glement is a necessary ingredient for excitation energy trans- 
port. 
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